function [rmin,rmax] = numericalBoundsPar(restr,irs,phi,opt,optimOptions)
% Function uses a constrained optimisation routine to find the upper and
% lower bounds of the identified set of the IRFs. This function requires
% the Parallel Computing Toolbox.
% Inputs are:
% - restr: structure representing restrictions
% - irs: structure containing IRs and possibly long-run cumulative IRs
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options
% - optimOptions: structure containing optimisation option

F = restr.F;
f = restr.f;
S = restr.S;
s = restr.s;
vma = irs.vmaGK;
Sigmatr = phi.Sigmatr;
Sigmatrinv = phi.Sigmatrinv;
H = opt.H;
ivar = opt.ivar;
jshock = opt.jshock;
signNorm = opt.signNorm;

N = size(Sigmatrinv,1); % No. of variables in VAR
mZ = sum(f); % No. of equality restrictions
mS = sum(s); % No. of sign restrictions (excl. normalisation)

% Draw opt.Qdraws initial values of Q for optimisation.
Q0 = drawQ0(restr,phi,opt);

% Set equality constraint for optimisation (Aeq*Q(:) = beq).
Aeq = zeros(mZ,N^2);
beq = zeros(mZ,1);
restrCount = 0;

for ii = 1:N % For each column of Q

    for jj = 1:f(ii) % For each restriction on ith column of Q

        restrCount = restrCount + 1;
        Aeq(restrCount,(ii-1)*N+1:ii*N) = F(restrCount,:);

    end

end

% Set inequality constraints for optimisation (Aineq*Q(:) <= b).
Aineq = zeros(mS+N,N^2);
bineq = zeros(mS+N,1);
restrCount = 0;

for ii = 1:N

    for jj = 1:s(ii)

        restrCount = restrCount + 1;
        Aineq(restrCount,(ii-1)*N+1:ii*N) = -S(restrCount,:);

    end

end

for ii = 1:N

    if signNorm == 0 % Normalisation on A0

        Aineq(mS+ii,(ii-1)*N+1:ii*N) = -Sigmatrinv(:,ii)';

    else % Normalisation on A0^(-1)

        Aineq(mS+ii,(ii-1)*N+ii) = -Sigmatr(ii,:);

    end

end

rmin = zeros(H+1,length(opt.ivar),size(Q0,3));
rmax = rmin;
LB = [];
UB = [];

parfor hh = 1:H+1 % For each horizon
    
    [rmin(hh,:,:),rmax(hh,:,:)] = genBounds(Q0,jshock,ivar,vma(:,:,hh),N,...
        Aineq,bineq,Aeq,beq,LB,UB,optimOptions)

end

rmax = -rmax;

% Take optima over different initial values.
rmin = min(rmin,[],3);
rmax = max(rmax,[],3);

end